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Abstract 

A first step toward a universal nuclear energy density functional based on low- 
momentum interactions is taken using the density matrix expansion (DME) of 
Negele and Vautherin. The DME is adapted for non-local momentum-space po- 
tentials and generalized to include local three-body interactions. Different prescrip- 
tions for the three-body DME are compared. Exploratory results are given at the 
Hartree-Fock level, along with a roadmap for systematic improvements within an 
effective action framework for Kohn-Sham density functional theory. 



1 Introduction 



Calculating the properties of atomic nuclei from microscopic internucleon in- 
teractions is one of the most challenging and enduring problems of nuclear 
physics. However, recent developments in few- and many-body physics to- 
gether with advances in computational technology give hope that controlled 
calculations of medium and heavy nuclei starting from a microscopic nuclear 
Hamiltonian will be forthcoming (see, for example, [T][2][5] ). Density functional 
theory (DFT), which is a self-consistent framework that goes beyond conven- 
tional mean-field approaches, offers particular promise for medium to heavy 
nuclei. The central object in DFT is an energy functional of the nuclear den- 
sities that would apply to all the nuclides. Phenomenological functionals have 
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had many successes but lack a microscopic foundation and theoretical con- 
trol of errors, such that extrapolations to the limits of nuclear binding are 
uncontrolled. 



Recent progress in evolving chiral effective field theory (EFT) interactions to 
lower momentum using renormalization group (RG) methods [¥|f5|l6|f7|l8|l9lll0|ll|r 
(see also p^|T^ ) makes feasible a microscopic calculation of a universal nuclear 
energy density functional (UNEDF) \l5l. The evolution weakens or largely 
eliminates sources of non-perturbative behavior in the two-nucleon sector such 
as strong short-range repulsion and the tensor force from iterated pion ex- 
change [9], and the consistent three-nucleon interaction is perturbative at lower 
cutoffs [7]. When applied to nuclear matter, many-body perturbation theory 
for the energy appears convergent (at least in the particle-particle channel), 
with calculations that include most of the second-order contributions exhibit- 
ing saturation in nuclear matter and showing relatively weak dependence on 
the cutoff [8j . These features are favorable ingredients for a microscopic Kohn- 
Sham DFT treatment [T6|17II18] . Indeed, Hartree-Fock is a reasonable (if not 
fully quantitative) starting point, which suggests that the theoretical develop- 
ments and phenomenological successes of DFT for Coulomb interactions may 
be applicable to the nuclear case for low-momentum interactions. 



A formal constructive framework for Kohn-Sham DFT based on effective ac- 
tions of composite operators can be carried out using the inversion method 
P^IEn]mi2^]|2Sll21l25]l25] . This is an organization of the many-body problem 
that is based on calculating the response of a finite system to external, static 
sources rather than seeking the many-body wave function. It requires a tractable 
expansion (such as an EFT momentum expansion or many-body perturbation 
theory) that is controllable in the presence of inhomogeneous sources, which 
act as single-particle potentials. This is problematic for conventional internu- 
cleon interactions, for which the single-particle potential needs to be tuned to 
enhance the convergence of the hole-line expansion [27][25] . but is ideally suited 
for low-momentum interactions. Given an expansion, one can construct a free- 
energy functional in the presence of the sources and then Legendre transform 
order-by-order to the desired functional of the densities. However, these are 
complicated, non-local functionals and we require functional derivatives with 
respect to the densities, whose dependences are usually only implicit. While 
this is a feasible program, it will require significant development to extend 
existing phenomenological nuclear DFT computer codes. 

We seek a path that will be compatible in the short term with current nuclear 
DFT technology but testable and systematically improvable. In this regard, 
the phenomenological nuclear energy density functionals of the Skyrme form 
have the closest connection to low-momentum interactions. Modern Skyrme 
functionals have been applied over a very wide range of nuclei, with quantita- 
tive success in reproducing properties of nuclear ground states and low-lying 
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excitations [29|30f31j . Nevertheless, a significant reduction of the global and 
local errors is a major goal [32]. One strategy is to improve the functional 
itself; the form of the basic Skyrme functional in use is very restricted, con- 
sisting of a sum of local powers of various nuclear densities [e.g., see Eq. ([T])]. 
Fits to measured nuclear data have given to date only limited constraints on 
possible density and isospin dependences and on the form of the spin-orbit 
interaction. Even qualitative insight into these properties from realistic micro- 
scopic calculations should be beneficial in improving the effectiveness of the 
energy density functional. 

A theoretical connection of the Skyrme functional to free-space NN interac- 
tions was made long ago by Negele and Vautherin using the density matrix 
expansion (DME) ^33ii3^|l35] . but there have been few subsequent microscopic 
developments. The DME originated as an expansion of the Hartree-Fock en- 
ergy constructed using the nucleon-nucleon (NN) G matrix [331311, which was 
treated in a local (i.e., diagonal in coordinate representation) approximation. 
In this paper, we revisit the DME using non-local low-momentum interactions 
in momentum representation, for which G matrix summations are not needed 
because of the softening of the interaction. When applied to a Hartree-Fock 
energy functional, the DME yields an energy functional in the form of a gener- 
alized Skyrme functional that is compatible with existing codes, by replacing 
Skyrme coefficients with density- dependent functions. As in the original ap- 
plication, a key feature of the DME is that it is not a pure short-distance 
expansion but includes resummations that treat long-range pion interactions 
correctly in a uniform system. However, we caution that the Negele- Vautherin 
DME involves prescriptions for the resummations without a corresponding 
power counting to justify them. 

The idea of using soft, non-local potentials in an expansion starting with 
Hartree-Fock was explored in the late sixties and early seventies (see, for ex- 
ample, Refs. [5Bll37|38j ). However, soft potentials were generally abandoned 
because of their inability to saturate nuclear matter at the empirical den- 
sity and energy per particleQ They have been revived in the context of low- 
momentum potentials (often referred to as "Viowfc") derived by transforming 
modern realistic NN potentials. The key to their success is the recognition 
that three-body forces (and possibly four-body forces) cannot be neglected. 
With lowered cutoffs, the density dependence of the three-body contribution 
drives saturation which accounts for the apparent past failure in nuclear 
matter when only two-body contributions were included. 

The present work is a proof-of-principle demonstration with a roadmap for 
future developments. We note the following omissions and simplifications. 



Calculations using hard NN-only interactions also fail to reproduce empirical 
saturation properties. 
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• We restrict ourselves to isoscalar (A^ = Z) functionals. This is merely for 
simplicity; generalizations to the full isovector dependence will be presented 
in the near future. We also defer inclusion of spin-orbit and tensor terms, 
which will require extensions of the DME treatment of Negele and Vau- 
therin [39] . 

• We work to leading order in the perturbative many-body expansion (i.e., 
Hartree-Fock). An upgrade path to include second order and beyond is 
described in Section [61 

• The form for the three-body force is limited to that of chiral N^LO EFT. 
This is consistent with current approximations used with low-momentum 
potentials, but will need to be generalized to accommodate evolved three- 
body potentials. 

• Pairing is essential for the quantitative treatment of nuclei, particularly 
unstable nuclei. The DME functionals described here can be adapted to in- 
clude pairing as done in conventional Hartree-Fock-Bogliubov phenomenol- 
ogy. However, a unified treatment is feasible with low-momentum interac- 
tions [iOlliT] . 

• There are unresolved conceptual issues for applying DFT to a self-bound 
system [121l43f44j that we will not address here (but which must be dealt 
with eventually). In addition, projection is not considered. 

Recently, Kaiser and collaborators have applied the DME in momentum space 
to a perturbative chiral EFT expansion at finite density to derive a Skyrme- 
like energy functional for nuclei [^461147] . Their analytic expressions for long- 
range pion contributions can be effectively applied in our formalism to avoid 
slowly converging partial-wave summations. However, we defer to future work 
a detailed description of this application and also comparisons with their re- 
sults. 

The plan of the paper is as follows. In Section we present the features of 
density functional theory needed in our treatment and discuss how applying 
the DME will lead us to a generalized Skyrme-like energy functional. In Sec- 
tion [3l we review the Negele/Vautherin derivation of the DME for non-local 
(in coordinate space) two-body potentials and make a direct extension to mo- 
mentum space. The result is a set of simple formulas for the basic coefficient 
functions in terms of integrals over partial- wave matrix elements of the MowA; 
potential. In Section HI we extend the DME to include three-body forces, re- 
stricting ourselves to local potentials of the form used in chiral EFT at N^LO 
(which is the form used in current approximations to low-momentum NNN in- 
teractions). We consider two prescriptions for the three-body part. We present 
some tests of the DME and sample results in Section [5], highlighting the effects 
of non-locality, the relative size of NN and NNN contributions, and the impact 
of different prescriptions for the NNN DME expansion. We conclude with a 
summary and roadmap for future calculations in Section [Hi 
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2 Density Functional Theory 



In this section, we give overviews of the standard Skyrme functional and the 
ideas behind Kohn-Sham DFT for nuclei that we need to set up the energy 
density functional calculations using the DME. 



2.1 Skyrme Hartree-Fock Energy Density Functional 



In the conventional Skyrme Hartree-Fock (SHF) formalism, the energy is a 
functional of the density p, the kinetic density r, and the spin-orbit density 
J. For simplicity, we restrict the discussion to N = Z nuclei here, so these are 
isoscalar densities only. This functional is a single integral of a local energy 
density, which depends in a simple way on these densities, such as 



Eshf[p, ^> J] = / + \top' + ^'^ap'^" + ^(3ti + 5t2)pr 

+ ^(9ti - ^t,){ypf - ^Wopv ■ J + ^(ti - 12) J'l . (1) 

Expressions for the Skyrme functional including isovector and more general 
densities can be found in Ref. [19]. The densities p, r, and J are expressed as 
sums over single-particle orbitals </>/3(x): 



p(x)^^|0,(x)p, (2) 

P 

r(x) = E|V0,(x)r, (3) 

P 

J(x)^5:0^(x)HVxa)0^(x), (4) 

where the sums are over occupied states and the spin-isospin indices are im- 
plicit. (More generally, when pairing is included with a zero-range interaction, 
the sums are over all orbitals up to a cutoff, weighted by pairing occupation 
numbers. This complicates finding the self-consistent solution significantly but 
is not important for our discussion.) The parameters to^^s, W^o, and a deter- 
mine the functional and are obtained from numerical fits to experimental data. 

Varying the energy with respect to the wavefunctions with Lagrange multipli- 
ers Bp to ensure normalizatioiH] leads to a Schrodinger-type equation with a 



^ Unconstrained variation of the orbitals is the usual textbook formulation of 



5 



position-dependent mass term [SUpS] : 

'^2M^^ + Ui^) + IWoVp --^Vxct] 0^(x) = M^) > (5) 

where ^ 

f/(x) = hop + ^hp^ + ^(3ti + 5t2)r + ^(5t2 - 9ti) - ^WqV • J , (6) 



the effective mass M*(x) is 

1 1 



+ 



2M*(x) 2M L16 16 



3 5 ■ 



P(x) , (7) 



and the Wo term is a spin-orbit potential (see Ref. [51] for details). The po- 
tentials in Eqs. (jHIl^dZD and the orbitals from Eq. (jS]) are evaluated alternately 
until self-consistency (see Fig. [TU]) . 

As we will see below, the DME energy functional for N = Z will take the 
same local form as -Eshf, 

^dme[p, r, J] = J d^RSuMEipin), r(R), J(R)) , (8) 

where the energy density function £^dme is evaluated with the local densities 
at R. We follow the Negele/Vautherin notation for £^dme and write [33] 

W = ^ + Mp] + B[p]r + C[p]\Vp\' + ■ ■ ■ , (9) 

where A, B, C are functions of the isoscalar density p instead of the con- 
stant Skyrme parameters, and we have suppressed terms that go beyond the 
present limited discussion. (When N ^ these are functions of the isovector 
densities as well.) Equation ([9]) implies that the DME form will be a direct 
generalization of the Skyrme functionals. 



2.2 DFT from Effective Actions 



Microscopic DFT follows from calculating the response of a many-body sys- 
tem to external sources, as in Green's function methods, only with local, 
static sources that couple to densities rather than fundamental fields. (Time- 
dependent sources can be used for certain excited states.) It is profitable to 

Skyrme Hartree-Fock [l8]. But this does not hold beyond Hartree level for a gen- 
eral microscopic DFT treatment with finite-range potentials, for which there is an 
additional constraint to the orbital variation |18j . 
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think in terms of a thermodynamic formulation of DFT, which uses the effec- 
tive action formahsm [52] apphed to composite operators to construct energy 
density functionals P^2Uf22j . The basic plan is to consider the zero tempera- 
ture limit of the partition function Z for the (finite) system of interest in the 
presence of external sources coupled to various quantities of interest (such as 
the fermion density). We derive energy functionals of these quantities by Leg- 
endre transformations with respect to the sources |53]. These sources probe, 
in a variational sense, configurations near the ground state. 

An analogous system would be a lattice of interacting spins, to which we apply 
an external source in the form of a magnetic field H [52]. The Helmholtz free 
energy F\li\ is calculated as the energy in the presence of the magnetic field 
and we determine the magnetization by a derivative with respect to the field, 
M{H) = —dF[H]/dH. It is often useful to reverse the problem, and ask what 
external field produces a specified magnetization. This leads us to the Gibbs 
free energy G[M], which we obtain by inverting M{H) to find H{M) and 
performing a Legendre transform: 

G[M] = F[H] + H{M)M . (10) 

Because H = dG[M]/dM and H vanishes in the ground state, G is extremized 
in the ground state (and concavity tells us that it is a minimum). If H is 
an inhomogeneous source, the formalism is generalized by replacing partial 
derivatives by functional derivatives and performing a functional Legendre 
transform. 



To derive density functional theory, we follow the same procedure, but with 
sources that adjust density distributions rather than spins. (We can either 
introduce a chemical potential or only consider variations that preserve net 
particle number. We implicitly assume the latter here.) Consider first the 
simplest case of a single external source J(x) coupled to the density operator 
p[x) = ip'''{x)ip{x) in the partition function 

Z[J] = e-^W ~ Tr e-^(-^+-^^ ^ Jv[^^]V[^] e-I^''+'^'^^ , (11) 

for which we can construct a path integral representation with Lagrangian C 
[52\ . (Note: because our treatment is schematic, for convenience we neglect 
normalization factors and take the inverse temperature /3 and the volume Q 
equal to unity in the sequel.) The static density p(x) in the presence of J(x) 
is 

p(x) = (p(x)), = ^ , (12) 
which we invert to find J[p] and then Legendre transform from J to p: 

T[p] = -W[J] + J d^x J(x)p(x) , (13) 
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with 



J(x) 



(5p(x) (5p(x) 







(14) 



Pgs(x) 



For static p(x), r[p] is proportional to the conventional Hohenberg-Kohn en- 
ergy functional, which by Eq. f|T^ is extremized at the ground state density 
(and thermodynamic arguments establish that it is a minimum j21j)J ^1 



Pi 



We still need a way to carry out the inversion from p[J] to J[p]', a general 
approach is the inversion method of Fukuda et al. p^1[20] . The idea is to 
expand the relevant quantities in a hierarchy, labeled by a counting parameter 
A, 



W[J, \]=Wo[J] + \Wi[J] + \^W2[J] + ■■■ , (15) 
J[p, X]=Mp] + AJi[p] + AV2H + ■ ■ ■ , (16) 

r[p, A]=ro[p] + Ari[p] + x't^Ip] + ■■■ , (17) 

treating p as order unity (which is the same as requiring that there are no 
corrections to the zero-order density), and match order by order in A to deter- 
mine the Jj's and Fj's. Zeroth order is a noninteracting system with potential 
Jo(x): 



ro[p] = -Wo[Jo] + / t/^x Jo(x)p(x) 



and 



5Jo(x) 

Because p appears only at zeroth order, it is always specified from the non- 
interacting system according to Eq. f|T9|) : there are no corrections at higher 
order. This is the Kohn-Sham system with the same density as the fully in- 
teracting system. 

What we have done is to use the freedom to spht J into Jo and J — Jq, 
which is essentially the same as introducing a single-particle potential U and 
splitting the Hamiltonian according to H = {Hq + U) + {V — U). Typically U 
is chosen to accelerate (or even allow) convergence of a many-body expansion 
(e.g., the Bethe-Brueckner-Goldstone theory j271l54f 28] ) . For DFT, we choose 
it to ensure that the density is unchanged, order by order. Thus, we need the 
fiexibility in the many-body expansion to choose U without seriously degrading 
the convergence; such freedom is characteristic of low-momentum interactions. 
(Note: If there is a non-zero external potential, it is simply included with Jq.) 



^ A Minkowski-space formulation of the effective action with time-dependent 
sources leads naturally to an RPA-like generalization of DFT that can be used 
to calculate properties of collective excitations. 
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We diagonalize Wo[<^o] by introducing Kohn-Sham orbitals (pi and eigenvalues 

[-VV2m- Jo(x)]0. = 5.0, (20) 

so that 

p(x)=x:i0.(x)p. (21) 

i=l 

Then Wq is equal to the sum of e^'s. The orbitals and eigenvalues are used to 
construct the Kohn-Sham Green's functions, which are used as the propagator 
lines in calculations the VTjf'^o] diagrams. Finally, we find Jo for the ground 
state by truncating the chain at Ti^^^, 

j^^W,^T,^J,^W,^T,^...^ W,^_ ^ r,_, (22) 

and completing the self-consistency loop: 

Calculating the successive Fj's, whose sum is directly proportional to the de- 
sired energy functional, is described in Refs. PU|21II55|I25] . 

When transforming from Wi to Fj, there are additional diagrams that take 
into account the adjustment of the source to maintain the same density and 
also so-called anomalous diagrams (these are two-particle reducible) . A general 
discussion and Feynman rules for these diagrams are given in Refs. [20ll21|23j . 
These two types of contribution cancel up through N^LO in an EFT expansion 
with short-range forces using dimensional regular izat ion [23], just as they do in 
the inversion method used long ago by Kohn, Luttinger, and Ward [SBjlFT] to 
show the relationship of zero-temperature diagrammatic calculations to ones 
using the finite-temperature Matsubara formalism in the zero-temperature 
limit. In the present application of the DME approximation to the effective 
action DFT formalism, they also cancel and so are omitted entirely. 

Note that even though solving for Kohn-Sham orbitals makes the approach 
look like a mean-field Hartree calculation, the approximation to the energy 
and density is only in the truncation of Eq. fl23|) . It is a mean-field formal- 
ism in the sense of a conventional loop expansion, which is nonperturbative 
only in the background field while including further correlations perturbatively 
order-by-order in loops. The special feature of DFT is that the saddlepoint 
evaluation applies the condition that there are no corrections to the density. 
We emphasize that this is not ordinarily an appropriate expansion for inter- 
nucleon interactions; it is the special features of low-momentum interactions 
that make them suitable. 

To generalize the energy functional to accommodate additional densities such 
as r and J, we simply introduce an additional source coupled to each density. 
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JnfR) 




+ 



Fig. 1. Schematic representation of Eq. ([27|) for a local potential, where the double- 
line symbol denotes the {6p/SJo)~^ term. 

Thus, to generate a DFT functional of the kinetic-energy density as well as 
the density, add ^^(x) 'Vip'^ • V?/' to the Lagrangian and Legendre transform to 
an effective action of p and r 



(24) 



r[p,r] = W[J,7]] - J d^x J{x)p{^) ~ J d^x7]{^)T{^) . 
The inversion method results in two Kohn-Sham potentials. 



^o(x) 



5rint[p, r] 



5p(x) 



and ?7o(x) 



(JFint [p,r] 



5tU) 



(25) 



where Tint = F — Fq. The Kohn-Sham equation is now [25 

1 



-V 



M*('x' 



■V - Jo(x) (pi = ei(f). 



(26) 



with an effective mass l/2M*(x) = 1/2M — i]q{x), just like in Skyrme HF. 
Generalizing to the spin-orbit or other densities (including pairing [ID]) pro- 
ceeds analogously. We note that the variational principle implies that adding 
sources will always improve the effectiveness of the energy functional. 

The Feynman diagrams for Wi will in general include multiple vertex points 
over which to integrate. Further, the dependence on the densities will not be 
explicit except when we have Hartree terms with a local potential (that is, 
a potential diagonal in coordinate representation). One way to proceed is to 
calculate the Kohn-Sham potentials using a functional chain rule, e.g.. 



MR) = 



int 



[P] _ 



6p{R) 



JMy) 5Jo(y) 



(27) 



and steepest descent [21] • This is illustrated schematically for a local interac- 
tion in Fig. [H We see that the Kohn-Sham potential is always just a function 
of R but that the functional is very non-local. If zero-range interactions are 
used, these diagrams collapse into an expression for Jo(R) that has no in- 
ternal vertices, but this is no longer true for diagrams with more than one 
interaction. Orbital-based methods take the chain rule in Eq. (1271) one step 
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further, adding a functional derivative of the sources with respect to the 0i's 
(and £i's); see Refs. [T8|58II59||60] for background on these calculations applied 
to electronic systems. Eventually, we plan to carry out such calculations to 
construct the full energy density functional. 

An alternative in the short term is to approximate Wi^t so that the dependence 
on the densities (rather than the sources or the orbitals) is explicit. This has 
two effects: the construction of the Fj from the Wi does not have additional 
terms and the necessary functional derivatives are immediate. An example of 
such an approach is the local density approximation (LDA). Here we go be- 
yond the LDA with the density matrix expansion (DME). By expanding the 
Wi about a "center-of-mass" R, we generate a local energy density that is a 
function of densities (p, r, . . . ) at R. We choose sources to match these den- 
sities and carry out the Legendre transformation implicitly; the end result at 
leading order is calculating Wi using density matrices built from Kohn-Sham 
orbitals. We are able to vary with respect to the orbitals because the con- 
straint of a multiplicative Kohn-Sham potential is built in. Then the resulting 
Kohn-Sham DFT has precisely the form of the Skyrme Hartree-Fock energy 
functional and single-particle equations. 



2.3 Low- Momentum Potentials 

The original DME application was based on a Hartree-Fock energy functional 
calculated with a G matrix, following the Brueckner-Bethe-Goldstone (BBG) 
method P71l54f28] . The latter involves infinite resummations of diagrams for 
nuclear many-body theory, as needed to deal with strongly repulsive poten- 
tials. In BBG there are two general resummations: the ladder diagrams into 
a G matrix and the hole-line expansion using the G matrix. Furthermore, to 
accelerate convergence of the hole-line expansion one needs to carefully choose 
a single-particle potential. This is problematic for the success of a Kohn-Sham 
DFT construction, for which the background field (which acts as a single- 
particle potential) has a separate constraint, namely to maintain the fermion 
density distribution. 

Renormalization group (RG) methods can be used to evolve realistic nucleon- 
nucleon potentials (e.g., chiral EFT potentials at N^LO), which typically have 
strong coupling between high and low momentum (i.e., off-diagonal matrix ele- 
ments of the potential in momentum representation are substantial), to derive 
low-momentum potentials in which high and low momentum parts are decou- 
pled. This can be accomplished by lowering a momentum cutoff A [H[5|6|7] 
or performing a series of unitary transformations that drive the hamiltonian 
toward the diagonal p^lllll2] . The UCOM transformations of Ref . [13] is an 
alternative to explicit RG methods. In all cases, we have a potential for which 
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only low momenta contribute to low-energy nuclear observables, such as the 
binding energy of nuclei. For convenience, we'll refer to any of these as Viowfc- 

We stress that evolving Viow k does not lose relevant information for low-energy 
physics, which includes nuclear ground states and low-lying excitations, as long 
as the leading many-body interactions are kept [TT]. The long-range physics, 
which is from pion exchange (and Coulomb), is preserved and remains local, 
while relevant short-range physics is encoded in the low-momentum potential 
through the RG evolution. Most important, for any Viowfc potential the obsta- 
cles from strongly repulsive potentials are removed. Hartree-Fock (including 
three-body interactions) saturates nuclear matter and G matrix resummations 
are not required (but may still be advantageous). Thus, we have a hierarchy 
suitable for DFT based on many-body perturbation theory. [Note: While the 
need for particle- hole resummations remains to be investigated for l^owfc po- 
tentials, results from the analogous UCOM potentials indicate perturbative 
particle- hole contributions for the energy [H].] 

While the evolution of Vfowfc potentials does not disturb the locality of initial 
long-range potentials, the short-range part becomes increasingly non-local. 
That is, in coordinate representation (r|V|r') has an increasing range in |r— r'|. 
Thus we must test that the DME is a good expansion for such non-localities. 

The interactions must include three-body (and higher-body) potentials, which 
should be consistently evolved with the two-body potential. These are not yet 
available (although SRG methods show promise of providing them in the near 
future [TI]|ll|12j ). and are instead approximated by adjusted chiral N^LO 
three-body potentials The vahdity of this approximation relies on the RG 
methods modifying only the short-distance part of the potential and is sup- 
ported by the observation that the EFT hierarchy of many-body forces appears 
to be preserved by the RG running [7]. The N^LO three-body potentials are 
local and we restrict our present investigation for now to this option. Given 
this microscopic NN and NNN input, we apply the density matrix expansion 
to derive an energy density functional of the Skyrme form. 



3 DME for Two-Body Potentials in Momentum Space 

In this section we derive the density matrix expansion for a microscopic DFT 
starting from low-momentum (and non-local) two-body potentials. From Sec- 
tion 12. 2^ the relevant object we need to expand is PFint, which is expressed 
in terms of the Kohn-Sham orbitals and eigenvalues that comprise the Kohn- 
Sham single-particle propagators. For Hartree-Fock contributions of the form 
in Fig. El^a), however, only the orbitals enter because the Kohn-Sham Green's 
function reduces to the density matrix. Similarly, higher-order contributions 
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such as the ladder diagrams in the particle-particle (pp) channel can also 
be put approximately into this form by averaging over the state dependence 
arising from the intermediate-state energy denominators. Therefore, while the 
results in this section are derived for the Hartree-Fock contributions to the 
functional, they can easily be generalized to include higher-order ladder con- 
tributions; this will be explored in a future publication. 

In essence, the DME maps the orbital-dependent expressions for contributions 
to Wint of the type in Fig. [2](a) into a quasi- local form, with explicit dependence 
on the local densities p(R), t(R), V^p(R), and so on. This greatly simplifies 
the determination of the Kohn-Sham potential because the functional deriva- 
tives of Pint can be evaluated directly. 

3. 1 Expression for PVhf 

Before presenting the details of the DME derivation and its application to non- 
local low-momentum interactions, it is useful to first derive in some detail the 
starting expression for H^hf, the Hartree-Fock contribution to Wint- This will 
serve to introduce our basic notation and to highlight the differences between 
most existing DME studies, which are formulated with local interactions and 
in coordinate space throughout, and the current approach, which is formulated 
in momentum space and geared towards non-local potentials. 

For a local potential, the distinction between the direct (Hartree) and exchange 
(Fock) contributions is significant, and is reflected in the conventional decom- 
position of the DFT energy functional for Coulomb systems, which separates 
out the Hartree piece. For a non-local potential, the distinction is blurred be- 
cause the Hartree contribution now involves the density matrix (as opposed 
to the density) and it is not useful to make this separation when the range of 
the interaction is comparable to the non-locality] ^ I Consequently, throughout 
this section we work instead with an antisymmetrized interaction. 

For a general (i.e., non-local) free-space two-body potential V , PVhf is defined 
in terms of Kohn-Sham states [Eq. (!20l) ] labeled by i and j, 

W^HF = ^E(^j|t>(l -Pi2)Kj) = ^E(u|V|zj) . (28) 

The summation is over the occupied states and the antisymmetrized inter- 
action V = V{1 — P12) has been introduced, with the exchange operator 

However, it is useful to separate out the long-distance part of the potential, which 
is local, and treat its direct (Hartree) contribution exactly. 
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Pi2 equal to the product of operators for spin, isospin, and space exchange, 
Pi2 = PcrPrPr- Notc that the dependence of VFhf on the Kohn-Sham potential 
has been suppressed. By making repeated use of the completeness relation 

11 = ^ f drlrar) {rarl , (29) 
VFhf can be written in terms of the coordinate space Kohn-Sham orbitals as 



ij {o-t} 

X 0*(ricriri)0j(r3cr3r3)0*(r2(T2r2)0j(r4cr4r4) . (30) 



From the definition of the Kohn-Sham density matrix, 

A 

p{rsa3T3, rio-iTi) = J2 (f>i i^i^^ in) 3(^3^3) , (31) 

i 

SO Eq. fl30|) can be written as 



^HF = 2 ^ jdri--- j dYi (rio-irir2cr2T-2|V|r30-3r3r4cr4r4) 

{o-r} 

Xp(r3(j3r3, ri(Tiri)p(r4cr4r4, r2cr2r4) 
= ]p:i{Ii2jdv,--- jdv, (rir2|Vi^>3r4)p^'nr3,ri)p('nr4,r2) , (32) 

where a matrix notation is used in the second equation and the traces denote 
summations over the spin and isospin indices for "particle 1" and "particle 
2" . Hereafter we drop the superscripts on V and p that indicate which space 
they act in as it will be clear from the context. 

Expanding the p matrices on Pauli spin and isospin matrices we have 
p(ri, Y2) = ^[po(ri, r2) + pi(ri, r2)r^ + So(ri, ra) • + Si(ri, ra) ■ arj , (33) 

where we have assumed the absence of charge-mixing in the single-particle 
states. The usual scalar-isoscalar, scalar-isovector, vector-isoscalar, and vector- 
isovector components are obtained by taking the relevant traces. 
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Fig. 2. (a) Schematic diagram for approximations to Wint that can be expanded using 
the DME. (b) Coordinates appropriate for the DME apphed to the Hartree-Fock 
potential energy with a non-local potential. 

A 

po(ri,r2) = Tr,.[p(ri,r2)] =5:0|(r2)0.(ri) , (34) 

i 

A 

pi(ri,r2) = Tr^r[/o(ri,r2)r^] = ^ 0|(r2)r^0i(ri) , (35) 

i 
A 

5o(ri,r2) = Tr..[p(ri,r2)a] = J] 0l(r2)a0,,(ri) , (36) 

i 

A 

Si{ri,r2) = TTar[p{ri,r2)aT^] = ^ 0j(r2)CTr^0i(ri) , (37) 

i 

where 0j(r) denotes a spinor with components 0j(rcrr). 

In this initial work we will only consider terms in the energy functional arising 
from products of the scalar-isoscalar (po) density matrices in Eq. fl32|) . which 
are the relevant terms for spin-saturated systems with N = Z. Thus, we will 
drop the subscript "0" on the density matrices from now on. 

After switching to relative/center-of-mass (COM) coordinates (see Fig. [2]) and 
noting that the free-space two-nucleon potential is diagonal in the COM co- 
ordinate, the starting point for our DME of the two-body Hartree-Fock con- 
tribution from a non-local interaction is 



32 



l^JdRdr dr'p{K + ^, R + ^)p(R " ^> R " ^)Tr.r [(r| V|r')] , 

(38) 



where V denotes the antisymmetrized interaction and the trace is defined as 
Tr,,[(r|V|r')] = ^ (^1X1^2X2! V(l - Pi2)|rViri(T2r2) . (39) 

{ar} 

The DME derivation of Negele and Vautherin (NV) [33] focuses on applica- 
tions to local potentials, which satisfy {r\V\r') = S{r — r'){T\V\r') . While the 
original NV work included coordinate-space formulas applicable for non-local 
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interactionq_l|, for low-momentum potentials it is convenient to revisit and 
extend the original derivation to a momentum-space formulation. We note 
that Kaiser et al. have shown how to use medium-insertions in momentum 
space in their application of the DME to chiral perturbation theory at finite 
density pIHGlliT] . 

For the momentum space formulation, we first rewrite the density matrices 
appearing in Eq. (138|) as 

p(R ± r72, R ± r/2) = p(R^ ± A/2, R^ t A/2) , (40) 
where the vectors appearing on the right-hand side are defined by (see Fig. [2]) 

R± = R±il], S = l(r' + r), A = l(r'-r). (41) 

Introducing the Fourier transform of V in the momentum transfers conjugate 
to S and A, 

q = k-k', p = k + k', (42) 
(where k', k correspond to relative momenta) gives 

= ^ / "^^z 9^ ' ^^3) 

where we have defined 

F(R, q, p) = J rfS dA e'"^-^ e'^'^ p(R+ - A/2, R+ + A/2) (44) 
xp(R- + A/2,R-- A/2) , 

and 

V(q, p) = 8 JdSdA e-'"^-^ e'^P'^ (S - A| V|S + A) . (45) 

The momenta q and p correspond to the momentum transfers for a local 
interaction in the direct and exchange channels. That is, the direct matrix 
element is a function of q and the exchange is a function of p. In contrast, for 
a non-local interaction the direct and exchange matrix elements depend on 
both q and p. This is the reason why we do not attempt to separate out the 
Hartree (direct) and Fock (exchange) contributions to Wuf, as is commonly 
done for local interactions. 



^ However, note that the final formulas for non-local potentials in Ref. [33] have 
numerous errors, which were not among those corrected in Ref. [34j . 
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The trace of Eq. (H5l) can be written in a more convenient form for our purposes 
as a sum over partial wave matrix elements, 

Tr.,[V(q, p)] = Svr ^ '(2j + l)(2t + 1) Pi{k ■ k'){klsjt\V\k'lsjt) , (46) 

Isj 

where the primed summation means that it is restricted to values where l + s+t 
is odd, with k = |(p + q) and k' = |(p — q). For simphcity we have assumed a 
charge-independent two-nucleon interaction, although charge-dependence can 
easily be included. 

3.2 Density Matrix Expansion 

The expression Eq. (1431) for VFhf is written in terms of off-diagonal den- 
sity matrices constructed from the Kohn-Sham orbitals. Consequently, the 
corresponding Enpfp] is an implicit functional of the density. The orbital- 
dependent Fhf requires the use of the functional derivative chain rule to eval- 
uate Ji(R) = 6Tnp[p]/6p(Il) in the self-consistent determination of the Kohn- 
Sham potential, which presents computational challenges and would require 
substantial enhancements to existing Skyrme HFB codes. 

Alternatively, we can apply Negele and Vautherin's DME to Wuf, resulting 
in an expression as in Eq. ([9]) with explicit dependence on the local quantities 
p(R),r(R), and|Vp(R)p, 

WnF = JdRiA[p]+B[p]r + C[p]iVpf + ---). (47) 

The starting point of the DME is the formal identity [33] 

p(R + s/2, R - s/2) = ^ 0*(R + s/2)0(R - s/2) 

a 

= [e-(^-^^)/^E</>*(Ri)'/>(R2)]^^^^^^^ , (48) 

a 

where Vi and V2 act on Ri and Ri, respectively, and the result is evaluated at 
Ri = Ri = R. We assume here that time-reversed orbitals are filled pairwise, 
so that the linear term of the exponential expansion vanishes. Hence, through 
second-order gradient terms the angular integral of the density matrix squared 
is equivalent to the integral of the square of the angle-averaged density matrix. 
In this way, the leading off-diagonal behavior of the density matrices in VFhf 
is captured by simpler expressions. 

The angle-averaged density matrix takes the form 
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p(R + s/2, R - s/2) = - / d cos exp s • ( Vi - V2) /2 p(Ri, R2) 



sinh[is 


Vi- V2I] 




V 


1- V2 





p(Ri,R2; 



(49) 



Rl=R2=R 



with s = |s|. Using a Bessel- function expansion (which is simply the usual 
plane- wave expansion with real arguments), 



1 1 °° 

xy X 



(50) 



k=0 



where Q is related to the usual Legendre polynomial by Q(-2^) = P2k+i{iz) / (iz) 
we can express the angle-averaged density matrix as 



p(R + s/2,R-s/2) 



J2{in + 3)j2n+l{skp{R)) 
n=0 

2\ " 



X Qn 



Vi - V2 
, 2A;f(R) 



p(Ri,R2) 



(51) 



where an arbitrary momentum scale A;f(R) has been introduced. Equation (!5T!) 
is independent of kp if all terms are kept, but any truncation will give results 
depending on the particular choice for kp. In this initial study, we employ the 
standard LDA choice of Negele and Vautherin: 



kp(K) 



(37rV(R)/2) 



1/3 



(52) 



Alternative choices for A;f(R) to optimize the convergence of truncated expan- 
sions of Eq. fl^T]) and to establish a power counting will be explored in a future 
paper. 

Following Negele and Vautherin, Eq. (!5T]) is truncated to terms with n ^ 1, 
which yields the fundamental equation of the DME, 



p(R + ^, R - ^) ^ Psl(/^f(R)s) p(R) (53) 

+s'g{kp{K)s)\]v'p{R) - r(R) + hp{Kfp{R) 
1-4 5 

where 

PsL(a;) = 3ji{x)/x , g{x) = 35j3(x)/2x^ , (54) 

and the kinetic energy density is t(R) = J2i |V0j(R)p. If a short-range in- 
teraction is folded with the density matrix, then a truncated Taylor series 
expansion of Eq. (1531) in powers of s would be justified and would produce 
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a quasi-local functional. But the local kp in the interior of a nucleus is typi- 
cally greater than the pion mass m.^, so such an expansion would give a poor 
representation of the physics of the long-range pion exchange interaction. 

Instead, the DME is constructed as an expansion about the exact nuclear mat- 
ter density matrix. Thus, Eq. fl53l) has the important feature that it reduces to 
the density matrix in the homogenous nuclear matter limit, pnm(R' + s/2, R — 
s/2) = pshikps) p. As a result, the resummed expansion in Eq. fl53l) does not 
distort the finite range physics, as the long-range one-pion-exchange contri- 
bution to nuclear matter is exactly reproduced and the finite-range physics is 
encoded as non-trivial (e.g., non- monomial) density dependence in the result- 
ing functional. The small parameters justifying this expansion emerge in the 
functionals as integrals over the inhomogeneities of the density. (See Ref. [2l] 
for examples of estimated contributions to a functional for a model problem.) 

In the case of a local interaction, the Fock term is schematically given by Wp ~ 
/ (iR(isp^(R+s/2, R— s/2)\/(s), so a single apphcation of Eq. fl53l) is sufficient 
to cast PVhf into the desired form. For a non-local interaction the calculation is 
more involved as two applications of the DME are required. Following Negele 
and Vautherin, we first rewrite the density matrices appearing in Eq. (1551) as 

p(R ± r'/2, R ± r/2) = p(R± ± A/2, R± t A/2) , (55) 
where the vectors appearing on the right-hand side are defined by (see Fig. [2]) 

R± = R±ls, S = i(r' + r), A ^ i(r' - r) . (56) 
To simplify the notation we define 

k^ = kp{R^), p±=p(R±), r±=r(R±), (57) 

and it is from now on understood that the functions without superscripts 
depend only on the center-of-mass vector R if the argument is not written 
exphcitly. 

The first application of the DME corresponds to an expansion in the non- 
locality A about the "shifted" COM coordinates R^, giving 

p(R + v'/2, R + r/2) = p(R+ - A/2, R+ + A/2) 

^ PsL{k^A)p+ + A'g{k+A) \]v'p+ - r+ + hfp+] . (58) 

L4 5 J 

Thus, we can expand the product of density matrices in Eq. (l38ll as 
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p(R + |, R + ^)p(R - I' - ^) = Psl(/c^ A)p+psl(A:p A)p- 

+ A2^(A;+A)psl(A:p A)p-[ivV " ^+ + P+ 

4 5 

+ A'gik^A)ps^ik^A)p-'[^\/'p- - r- + ^fcp'p-] , (59) 
where we have dropped terms quadratic in the gradient. We then define 



«(p^)=PsL(fcM)P^, (60) 

and use Eq. fl50l) to perform a second density matrix expansion on a{p'^)a{p~) 
in S about R, 

5]2 g 
a(p+)a(p") ^ Psl(A;fS) + — ^(A;FS)[aV^a - |Va|^ + -Ajfa^] . (61) 

^ 



From a Taylor expansion of psl(^fS) and glkpTj) it is evident that the (/cpS)^ 
coefficients of exactly cancel each other. Because we desire a final expres- 
sion that reproduces the exact nuclear matter limit (and the presence of the 
PshikpT,) term spoils this limit), we follow the philosophy of Negele and Vau- 
therin and use this leading cancellation to motivate a different rearrangement 
and truncation of Eq. flSTj) such that 



5]2 

a(p+)a(p-) + — ^(A;FS)[aV^a - |Va|2] . (62) 

The freedom to rearrange the expansion as in the last equation stems from 
the fact that the restriction of Eq. ( 15T1) to n ^ 1 terms gives a truncated ex- 
pansion in powers of S^. The neglected terms, starting with S^, involve higher 
derivatives of the density. But having neglected these terms, retaining the 
other (and higher) contributions that are summed in g{k-pT?) is somewhat 
arbitrary. Therefore, Negele and Vautherin argue that it is advantageous to 
use this arbitrariness to "reverse engineer" the expansion so that the exact 
nuclear matter limit is always exactly reproduced by the leading term [33] . 
We emphasize that this is a prescription without established power count- 
ing or error estimates, which must be assessed in future work. As we show 
in Section [5l different prescriptions can lead to significant changes in nuclear 
observables. 

The gradient terms in the above equation can be evaluated with the aid of the 
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chain rula ^ I 

V«(p) = Vp^ , VMp) = VV|^ + |Vpr0 . (63) 

RecaUing that we define the local Fermi momentum as kp = (Svr^p)^/'^, we can 
explicitly evaluate the first and second derivatives of a, 

_=,„(fcrA), _ = -^,,(*,A). (64) 

Pulling it all together, the product of density matrices in Eq. (|38|1 are approx- 
imately given in terms of local quantities by 



p(R + ^, R + ^)P(R - - ^) ^ pkikFA)p' + ^J^'gik^J:) 

X (pVVpsl(A:fA)jo(A;fA) - \V pf[foikFA) + jf{kpA)]] 

+ 2A2^(fcFA)psL(A:FA)(ipVV - pr + ^A:2p2) . (65) 



3.3 Evaluation o/F(R, q, p) and the DME coupling functions 



In the momentum space expression for Whf, it remains to evaluate the Fourier 
transforms defined in Eq. (jHj) for the expanded density matrices in Eq. (1651) . 
Identifying the terms in Eq. (H71) that give the DME functionals A[p], -B[p], 
and C[p], we have 



F(R,q,p) 
F(R,q,p) 
F(R,q,p) 



A-K r 6 

^ = -(27r)35(q)^/2(p)pr, 



(66) 
(67) 



27r 



/3(g)/5(p)|Vp|^+ (27r)^(5(q)-^/2(p) 



L5 



+ 



ft/ TP 



(6J 



where p = p/kp etc., and the R-dependence of kp, p, and r has been sup- 
pressed. The functions Ij{p) and Ij{q) are simple polynomials (and theta func- 
tions) in the scaled momenta p and q: 



^ Note that the equations here assume the canonical choice of kp = (37r^/3/2)^/^. 
Alternative choices for kp, such as the one proposed by Campi and Bouyssy [61j 
where ^f = /cf(P; V^p, r) will generate additional terms by the chain rule. 
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r 377- 

h{p) = Jx^dxUpx)pl^ix) = —il6-12p + p')ei2-p), (69) 

l2iP) = J x'^dxjoipx) psLix)g{x) 

= ip^- 18p^ + 40p2 - 24p) e{2 - p) , (70) 

128 

r 3^7r 
h{q) ^ / x' dxjoiqx) g{x) = - — {bf - 3) ^(1 - q) , (71) 



hip)= I x^dxjo{px)jo{x)psL{x) = ^i2-p)9i2-p) , (72) 



h{p)= I x'dxjo{p-x)[fo{x)+jKx)] = ^{A-p')d{2-p) . (73) 



Note that the trivial angular dependence of Eqs. fl69|) - fl73|) is a consequence 
of the angle averaging that is implicit with each application of the DME. 

With the aid of Eqs. (l66l) -( |73l) . we can now obtain explicit expressions for 
the A, B, and C coupling functions by grouping terms appropriately and 
performing the relevant angular integrals. The expressions for A and B follow 
immediately and are given by 



B[p] 



^ / V^dpTT,r[V[Q,p)]{h{p) + -h{P)) 

5 



IGvr/cp Jq 



2kp 



P^dpTl^r[V{0,p)]l2{p) 



(74) 
(75) 



where Tra-r[V(0, p)] is given by a simple sum of diagonal matrix elements in 
the different partial waves, 

Tr..[V(0,p)] =87r^ ' (2j + l){2t + 1) {^lsjt\V\^lsjt) . (76) 

Isj 

The primed sum is over all channels for which / + s + t is odd. 

The contributions to VThf that have gradients of the local density take the 
form 



HF 



|Vp|2 



dR (Cv2pVV(R) + qvpp|Vp(R)|' 



(77) 



We can perform a partial integration on the V^p terms to cast them into the 
canonical form proportional to only |Vpp; that is, 



HF 



|Vp|2 



dR\Vp{R)\' 



so that 



d 



(78) 



(79) 
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In practice it is efficient and accurate to calculate the derivative in Eq. (179!) 
numerically rather than analytically. 



The expressions for C|vp|2 and CyZp are obtained by substituting the relevant 
terms in F(R, q, p) [see Eqs. fl66l)- fl67j) ] into Eq. fH3|l and performing the 
angular integrals, 



= -^^rnJ ^'^W P'dph{q)hiP)Vav{q,p) , (81) 

lOTT rh-p Jo Jo 



2fe 



F 



p^rfp/2(p)Tr,,[V(0,p)] 



+ I6^i, ^''^'^i) P'^P^3(^)^4(p)Vav(g,p), (82) 

where Vav(5',p) is the angle- averaged interaction, 

Va.{q,p) = \j rf(cose)Tr..[V(q,p)] , (83) 

and V(q, p) is given by Eq. fl46p . Note that care must be taken in the eval- 
uation of dC^2p/dp if the vertex V(q, p) is density-dependent or if the local 
Fermi momentum is not taken to be fcp = (37r^p/2)^/^. 



4 DME for three-body potentials in momentum space 



In this section we extend the DME as applied to the Hartree-Fock energy 
to include three-body force contributions. The low-momentum interactions 
currently in use do not yet include consistently evolved three-body forces be- 
cause of technical difficulties in carrying out the momentum-space evolutioiTH. 
Therefore, as an approximation to the evolution, two short-distance low-energy 
constants in the leading chiral three-body force (this is N^LO according to the 
power counting of Refs. [63|6l] ) are fit at each cutoff to properties of the 

^ However, the recent application of similarity renormalization group (SRG) meth- 
ods to inter-nucleon potentials provide a computationally feasible path to the 
momentum-space evolution of many-body forces |10|62j . 
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triton and ^He to determine the three-body force. In the present work, we 
will use this force exclusively and postpone the treatment of general non-local 
three-body forces, as will be produced by an SRG evolution. 



4-1 Wup for local three-body forces 

The Hartree-Fock 3NF contribution to the total energy is given by 

Wg''^ = ^j:{tjk\VA^23\tjk) , (84) 

ijk 

where the summation is over the occupied Kohn-Sham states and the operator 
^123 is the (un-normalized) three-nucleon antisymmetrizer 

Al23 = (1 + P13P12 + P23Pl2)(l - P12) 
= (l + P13P23 + Pl2P23){l - P23) 

= (1 + P23^'l3 + Pl2Pl3)(l - ^^13) • (85) 

Decomposing the three-body potential in the standard fashion [65] , 

V = + + V^(3) ^ (gg^ 

where 1/^*^ is symmetric under j ^ k, we can write the full interaction in 
terms of one component 

V = + P23Pl3V^'^Pl3P23 + P23Pl2V^'^Pi2P23 , (87) 

and so on. This allows us to simplify Eq. flMl) by using 

V^A23 = (1 + ^^23^13 + P23Pl2)V^'Uu3 , (88) 

the cyclic nature of the trace along with (1 -|- P23P12 + -P23-Pi3)-4i23 = 3^123, 
and other permutation operator identities to obtain 

ijk 
1 ^ 

= 2 T.^^3k\V^^\l + ^^23^12 + P13P12 - P12 - P23 - Pi3)\tjk) 

ijk 
1 ^ 

= 2 Y.{^jk\V^'\^ + 2P23P12 - 2Pi2 - P23)\^3k) . (89) 

ijk 
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Because the leading chiral EFT 3NF has a vanishing direct piece, there are 
only three independent contributions to Wuf that need to be evaluated: one 
double-exchange term involving two permutation operators and two single- 
exchange contributions. Writing Eq. (1891) in terms of density matrices and 
separating out the scalar-isoscalar contributions to wj^-^^ arising from single- 
exchange terms gives 



P(X2, Xi)p(xi, X2)p(x3) Tri23[V^^^Hxi, X2, X3)P] 



(JTl 

12 J 



1 



+ 2^'^^3, X2)p(x2, X3)p(Xi) Tri23[V^^^^ (Xi, X2, X3)P27] 
~ J rfXi C/X2 rfXg |p(x2, Xi )p(xi , X2)p(x3) 

/ ^|^e-^--(---)e---(-x3)Tr,,3[v(i)(q„ q3)P^ 

+ ^P(x3,X2)p(x2,X3)p(xi) 

^^^e-^-(---)e-^^3-(x,-x3)Tr^^3[t7(i)(q,, q3)p-]) | 



(90) 



where Tri23 = Tra.iriTro-2r2Tra-3r3 and a local 3NF has been assumed. Similarly, 
the scalar-isoscalar contributions to VFhf arising from the double-exchanges 
are given by 



W^HF ^ = -^J C?X3 p(Xi, X2)p(x2, X3)p(x3, Xi) 

where the Fourier transformed 3NF components are defined by 

(kik2k3|y«|k;k;,k'3) = (^)'5(q, + q2 + q3)y«(q2,q3) . (92) 

Here Q is the volume (which drops out of all final expressions) and qj = kj — k- 
is the momentum transfer. 

As discussed above, we approximate the RG evolution of the 3N force with the 
leading-order chiral 3N force, which is comprised of a long-range 27r-exchange 
part Vc, an intermediate-range Ivr-exchange part Vd and a short-range contact 
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Fig. 3. The chiral three-body force at N^LO according to the power counting of 
Ref. [6^, which has a long-range 27r-exchange part Vc (left), an intermediate-range 
Ivr-exchange part Vd (middle), and a short-range contact interaction Ve (right). 



interaction Ve [631IM] , see Fig. [31 The 27r-exchange interaction is 



2/ J {ql + ml){q^^ + ml) ' ' ' 



(93) 



where is defined as 



pa/3 _ ra/3 



4cim^ 2c3 



C4 



+ 2. 7i e"''" ■ (q^ X q.) ' (94) 

7 J TT 



while the lyr-exchange and contact interactions are, respectively, 



qa cd o-j ■ qj 



CE , ^ 



Tj) {at ■ qj) , 



(95) 
(96) 



In applying Eqs. dnSD-dHSD, we use = 1-29, U = 92.4 MeV and = 
138.04 MeV and the Cj constants extracted by the Nijmegen group in a partial 
wave analysis with chiral 27r-exchange [66]. These are ci = — 0.76GeV~^, 
C3 = — 4.78GeV~^ and C4 = 3.96 GeV^^. Fit values for the cd and ce low- 
energy constants consistent with a sharply cutoff low-momentum potential 
are tabulated in Ref. [671 for = 700 MeV. 



From the previous general expressions for W^p and W^p, we need to evaluate 
the spin-isospin traces Trias [V'^^^Pfs^], ^^i23[V^^^ P23"'], and Trias [V'^^^PaT^i" 
For the single-exchanges we find 



12 J 
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Tri23[Fi'^(q2,q3)P27]=48 



Tri23[Vi'Hq2,q3)P, 



ICTTl 

23 J 



Tri23[K^'nq2,q3)P^ 



23 J 



-48 



4 



:48f^f q2-q3 



^2fJ {ql + ml){ql + ml 
X 



4cim2 2C3 
—r2 H ^ q2 ■ q3 

J TV J TV 



Tri23[V^i'Hq2, q3)Pr2l = Tri23[Vi'Hq2, qs)^^] 

= Tri23[v;^'Hq2,q3)AT] = o, 



(97) 
(98) 

(99) 
(100) 



while the various double-exchange terms give 



(TTl 
12 J 



Tri23[V#^(q2,q3)P27^i 

Tri23[^i'^(q2,q3)P27AY 
Tri23[V;«(q2,q3)P27^r 



12 



12 



-12 



5'A CD 



12 (^V- q2-q3 



(101) 
(102) 



X 



Aciml 2c3/ \ 

+ — (1 + C4/C3j q2 ■ q3 



/•2 /■2 
\of_J f2 



VJ PAql + mlM + ml) 



. (103) 



Note that for the Ve and Vd terms, it is not necessary to treat separately the 
single- and double-exchange contributions because their structure is identical 
due to the nature of the zero-range three- and two-body vertices. Substituting 
the spin-isospin-traced interactions into Eqs. fl90l) - fl9Tl) and simplifying gives 



-^9E J rfx [p(x)]3 , 

go J rfX2 rfxa [p(x2,X3)]V(x2) 



3 

16 



(104) 
(105) 



dq: 



{2tt 



3 ^-iq3-(x2-X3) ?3 



qi + ml 



where gE = ce/JAx and gn = {gA/^fl) {cn/pAx)- Similarly, the single- and 
double-exchange contributions from the 27r-exchange 3NF are given by 
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{lx,c) 



HF 



3 r 

--Qc J dXi rfX2 dXg p(Xi)p(x2, X3)p(x3, X2 

dq2 dqs 



X 



(27r)6 



g-iq2-(xi-X2)g-iq3-(xi-X3) 



J IT 



^2C3 

+ -7^ q2 • q3 

/ IT 



q2 ■ qs 



{q2+ml){qi + ml) 
(106) 



and 



Tj/(2a;,c) 
t'KHF 



3 

16 



gc J dxi dX2 dxs p(xi, X2)p(x2, X3)p(x3, Xi) 

/■dq2C 
""J f27r 



^q3 iq2.(xi-X2) „-iq3-(xi-X3) 

(27r)6 



q2 • qs 



(g| + m2)(g| + m2) 



4cim^ 2c3 , . , 

X 75 h ^(1 + C4/C3) q2 ■ q3 

J IT J TT 

2C4 



/2 (g| + m2)(g2 + ^2); ' 



(107) 



where Qc = {qaI'^^UY. 



4-2 D-term 



As with the nucleon-nucleon contributions to Whf, it is convenient to re- 
cast the 3NF Hartree-Fock expressions into momentum space. Changing to 
relative/center-of-mass coordinates (R = (x2 + X3)/2, r = X2 — X3), the Itt- 
exchange 3N Hartree-Fock contribution becomes 



W^F^^go J cZRdr[p(R + r/2,R-r/2)]V(R + r/2) 



d(l _-iq.r <f_ 



(27r)3 g2 + ^2 

where we have defined 

F(R, q) = Jdr e-^*»""[p(R + r/2, R - r/2)] V(R + r/2) . (109) 



Applying the DME separately to the product of non-local and local densities 
in F(R, q) yields 
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[p(R + r/2, R - r/2)f ^ p^i.{k^r)p + r''g{k^r) \\pV^p - 2pr + hip" 



(110) 



and 



rl 3 

p(R + r/2) ^ psL(A;Fr)p + r'^gik^r) -VV + 

1-4 5 



(111) 



Combining the two expansions and dropping terms of higher order in the 
DME, we find 



[p(R + r/2, R - r/2)] V(R + r/2) ^ pUkFr)p^ 

"3 6 
+ r^g{k^r)psi^{k^r) -^P^V^P - 2pV + -kjp^ 



(112) 



where the R-dependence of kp and the local densities has been suppressed. 
Evaluating the Fourier transform defined in Eq. f ll09p using the approximate 
DME expressions and grouping terms according to which coupling function 
contribute gives 



F(R,q)|^ = 4vr(i^j [h{q) + -h{q) 



3 , 



F(R,q) 



U5 



3tt 



F{n,q)^ = ^p'V'ph{q) 



(113) 
(114) 
(115) 



where the integrals Ii{q) and hiq) were defined in Eqs. ( I69|) -(170l) and q = 
q/kp. Together with Eq. fllOSp . we obtain the Ivr-exchange 3NF contributions 
to the EDF coupling functions 



Ad[p] 
BdIp] 
CdIp] 



^^\d fdq^^[h{q) + ^hm 
J q^ + 5 



Sukl 



3p2 



47r/c| 



9d 



'^I^TZ — i"^2(g) 



327r 



dp \/cp 



^^^T — ^^2(g) 



(116) 
(117) 
(118) 
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4-3 c-term single- exchange 



Starting from the single-exchange HF contribution of the 27r-exchange 3NF in 
Eq. (11061) . we first change to Jacobi coordinates, 

r23 = X2 - X3 , ri = X2 - ^ (X3 + Xi) , R- = ^ (xi + X2 + X3) , (119) 

followed by the change of momentum variables q = |(q2 — Qa) and p = q2 + q3- 
The result is 



--g^ J dKJ -^^-^Fi^(R,p,q)K;iC3(p,q) 



(120) 



where Fi3;(R, p, q) is the Fourier transform of the product of density matrices, 



Fi,(R, P, q) = / dr, dr,^ e-'^'^' e'^-'^' p(R + 2ri/3) 

X [p(R - ri/3 + r23/2, R - ri/3 - r23/2)]2 

and Vcj^cs (p, q) is defined as 



:i2ii 



Vcicaip, q) 



q2 ■ qa 



iciml 2c3 



(qi + m2)(q2 + m2; 
with q2 = p/2 + q and q3 = p/2 — q. 
Referring to Eq. (11211) . we first expand p(x2,X3) as 



f2 + 72"^2 ■ qs 



(122) 



p(x2, X3) = p(R - ri/3 + r23/2, R - ri/3 - r23/2) 

^psL(A;F(R")r23)p(R") + r^3^(fcF(R")r23) 



X 



t(K- 



:fci(R-)p(R- 



-4 ^ ' ' 5 

where R^ = R — ri/3. Performing a subsequent expansion about R gives 



(123) 



rl 3 
p(x2, X3) ^ PsL(/^Fr23) P + r23fi'(A;Fr23) [^^V " ^ + ^^fP 



+^r2^(A;Fr23)V='(psL(fcFr23)p) , (124) 



9 

where the second application of the DME has been modified slightly to ensure 
the leading term is exact in the nuclear matter limit. Similarly, the diagonal 
density p(xi) is expanded as 
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2 4 
p(R + -ri) ^ p + -rlgikFri)V^p . 



(125) 



Therefore, to second order in the DME we obtain 



p(R + 2ri/3) [p(R - ri/3 + r2s/2, R - ri/3 - r2s/2)f ^ [psLihr23)]^p' 

rl 3 

+ 2rls9{kFr23)pSLikFr23) [^P^^V - pV + -fcpp^ 

1 4 
+ -r2(7(A;Fri)V2(psL(A:Fr23)p) + ^rjg {kFn)pl^{kFr 23) p^V^ p ■ (126) 



For the usual LDA choice for fcF(R), the V^(pslp) term evaluates to 



V^(pSL(A;Fr23)) 



d 

{pSL{kFr23) + P^PsL(A;Fr23)}vV 

+ {2—pSL{kFr23) + P^PSL(/i;F?^23)}|Vpp 



dp 



dp^ 



(127) 



which suggests a grouping of terms in Eq. fll26p according to which coupling 
function they contribute to, 



p-p 



p-p 



p-p 



6 

PsL(^Fr23)p^ + ■^rlr^g{.kFr23)psi.{kFr23)klp^ , 

^ -2r235'(/i;F^23)PSL(/i;F?^23)P^1" , 

1 2 

2^235'(^Fr23)PSL(A;Fr23)p^ + g'^lfi'(A;Fri)psL(/i;Fr23)P^ 



d 

X {3psL(A;Fr23) +P^PSL(A;Fr23)} 

■2 



+ 



^ri5((fcFri)PSL(A;Fr23)P^ 



X {2 — P5(A;Fr23) +P^PSL(/^Fr23)} 



9p^ 



Vp 



(128) 
(129) 



(130) 



Evaluating the Fourier transform in Eq. (11211) gives 
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2p3 



6 



Fi^(R,p,q) 
Fi^(R,p,q) 



Fi.(R, P, q) ^ = (27r)^5(p)-^ + -h{q) 

ftp 

(2vr)^5(p)f^J2(g) 

Hp 

+ ^^^3(p)(/i(g)--/6(g) 



(131) 
(132) 



c 



327rV 
9H 



/3(p){/6(g) + ^/7(g)-k(g)} 



Vp ,(133) 



where /1-/5 have been defined in Eqs. (!69l) - (fr3l) and the new integrals Iq-Is 
are defined as 











hip)^ 






37r 




32^ 



37r 

x^c^x jo(pa;)psL(x)j2(a;) = — (8 - 8^ + ^^) 6'(2 - p) 
Sn , 



■ip 



-8 + 40p - 36p' + 5p*) e{2-p) . 



(134) 
(135) 

(136) 



With exphcit expressions for the DME approximation to Fi3;(R, p, q) in hand, 
all that remains is to insert Eqs. f ll31l) -( fT33|) into Eq. (11201) and group terms 
accordingly. The A [p] and B [p] coupling functions follow immediately and are 
given by 



(137) 

(138) 



The derivation of the C[p]2^ coupling is a bit more comphcated because we 
must first partially integrate all V^p terms. Writing the gradient contributions 
to as 



HF 



|Vp|2 



c^-.,vV(R) + qV%p|Vp(R)h 



dR 

V ^ 

dR|Vp(R)p 



dp 



V2p 



(139) 
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we obtain 



C[p]\l = CfV^p - . (140) 



Comparing to Eqs. fll20p and (11331) we find 



^^^P^" ^ ~ Un^Pkl I ^'"^^ ^'"^^ ^'^^^ 

x{-hiq)-^Ijiq) + \lsiq)} (141) 



and 



327r/2A;| 

2 2 1 

- ^^1^ / p'dpq'dqV,,,,ip,q)h{p){hiq) - -Ie{q)} , (142) 
where the angle-averaged interaction V^cics {p, <?) is defined as 

Vc^ca {P,q) = ^j d cos ^ Vc,c3 (p, q) • (143) 

4-4 c-term double exchange 



The double exchange contribution from the c-term is given in Eq. (19 ip . Since 
this involves a product of three off-diagonal density matrices, the DME is 
significantly more involved than for the other 3N contributions. In order to 
assess the sensitivity to the details of the (non-unique) DME prescription, 
we consider two different expansion schemes for these contributions, which 
we denote by DME I and DME II. We expect the differences between the two 
schemes should be "small" if the master formula Eq. (153!) is indeed a controlled 
expansion, and if results are insensitive to the different angle-averaging used 
in the two schemes. 



4.4.1 DME I 

We start by noting that repeated application of the master formula Eq. (153!) 
factorizes the three-body center-of-mass and relative coordinate dependence 
as 

p(xi, X,) = J2 Mrm, rij)Oi{R) . (144) 
I 

where Oi{R) is some monomial of the local densities and i,j,m are a permu- 
tation of 1, 2, and 3. The relative coordinate functions can be written in terms 
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of their Fourier transforms, e.g., 

Hkrr. = / ^|^e^'^™^™e^'^-'-'^A(r^, r,,) . (145) 

Expanding the appropriate set of Jacobi coordinates for each density matrix, 
Eq. ( 19T|) can therefore be written as 

^^'^ ''^ ^ (2ttW ^16 ^ / '^^^ d^3d(l2 dqs ^,^304 (p, q) 

V ' ijm 

X A,(A;i, h3)X,ik2, h,)XUh, ku)OiiR)0,iR)OkiR) 

^ g-j(ki ri+k2 r2+k3 r3+ki3 ri3+ki2 ri2+k23 r23+q2 ri2-q3 ri3) (145) 

with Qc = {gA/'^fn)'^ and where Pk denotes an integration over all variables 
of type km and kj^ and 



M ( \ 12-43 r 4:Ciml 2c3 / x 

?M (U7) 

Now choose one set of Jacobi coordinates, e.g., V2 and ris and rewrite Eq. (11461) 
in terms of these alone 



13 1 
ki • ri — > ki • (--r2 - -ris) , kaa • rss — ^ kss • (rg - -rig) , 

kg • rg > kg • ("^fs + ^Fig) , ki2 • > ki2 • ("Fs - ^Fig) , 

q2-ri2 — ^q2 • (-r2 - ^ria) . (148) 
We obtain as our final result 



xOi{R)Oj{R)Ok{R)V,,,,,,{lCi,lC2) , (149) 

with 



/Ci = ks - ^ki - ^kg + k2g - ki2 , 

/C2 = kig-iki + k3-k23. (150) 
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Now let us consider the particular form of the functions appearing in the 
integrals. We expand each density matrix as in Eq. fll24p and use Eq. ( 16^ to 
evaluate the V^{pslp) term: 



p(xi,X2) ?5ip 



Pshikpru) + rl^gi^kpr 12) -k], 

5 



+ r 



rl^g{kFr 12) 



+V^p^-^g{kFr^)]Q{kFri2) + '^-fg{kFri2] 



+ivpr 

This leads us to define 



{kFri2? 



PShikFTu) 



(151) 



Ai(r3,ri2) = (pshikFru) + rl^g {kFr 12) -k^ 
>^2ir3,ri2) = -rl29ikFri2) , 
A3(r3,ri2) = -^g{kFn)jo{kFri2) + -^gikFru) 



>^4(r3,ri2] 



r 



{kFruY 



PshikFru) ■ 



(152) 
(153) 
(154) 

(155) 



We obtain the A-term by inserting the relevant functions into Eq. 01491) 



16(27r; 
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mAi(A;3,fci2)Ai(fc2,fci3)Ai(A;i,A;23)KiC3C4(/Ci,/C2) , (156) 



with 



Al(fc3, ki2) 



Qtt^ 2l7r2 



U3 
hp 



Ik^F 



(34 - ^kl^) e{kF - ki2){27c)'5^^\k3) 



X,{k,2){27rf6^'\k, 



(157) 



Integrating over the (5-functions leads to 



16(27r)s 



dkudkis dk23 Ai(A;i2)Ai(A;i3)Ai(A;23) 



X ¥^,0304(^23 - ki2, ki3 - k 



23 j 



:i58) 



The B-term is proportional to r 
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^[p] = Y^^^ / ^^^=1^3^4(^15 ^2) ^A2(A;3,A;i2)Ai(A;2,A;i3)Ai(/ci,A;23) 

+Ai(A;3, ki2)X2{k2, ki3)Xi{ki, /C23) + Ai(A;3, ki2)Xi{k2, ki3)X2{ki, ^23)^ , 

(159) 

with 

A2(A;3, k,2) = -^(34 - 5kf2)eikF - k,2)i2nf6^'\ks) 

= A2(A;i2)(27r)35(3)(k3) . (160) 

Integrating out the 5-functions gives 



B[p\ = ^Q^2^Y j ^ki2 ciki3 rfk23V;,c3C4(k23 - ki2, ki3 - k23) 

X ^A2(/ci2)Ai(A;i3)Ai(fc23) + Ai(A;i2)A2(fci3)Ai(fc23) 

+Ai(A;i2)Ai(A;i3)A2(A;23) ) , (161) 



The calculation of the relevant integrals for the C-term is more involved. We 
first consider on the integral for the coefficient of |Vpp 



C|vpp = iQi^2'Ky^ J ^'^^=1^3^4(^15 ^2) (Ai(/i;3,A;i2)Ai(A;2,A;i3)A4(A;i,A;23) 

+Ai(A;i, k23)Xi{k2, ki3)X4{k3, k^) + Xi{ki, /c23)Ai(/i;3, ki2)X4^{k2, kia) 

(162) 

with 

X^iki, k23) = -^y;^(34 - 5kl)Q{kF - ki) 

X f r^^(fcF - A;23) + i^J-^ikp - k23)] . (163) 

\kFk23 k23dk23 J 

Let us focus our attention on the first term in Eq. (11621) . This term contains 
a factor 

(^5{kF - k23) + :^S{kF - k23)] 7^V;,e3C4 , (164) 
\Kf Cl/v23 / fc23 
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which simphfies to 



^ 5{kF-k2^)V2x-^5{kF-k2^):^V2x. (165) 



fc23fcF «;23 dfc23 

after partial integration. 

The second term has the coefficient V^p 



Cvv = y|§^ / ^'^^-i^3C4(/Ci,/C2)(Ai(A;3,A;i2)Ai(A;2,A;i3)A3(A;i,A;23) 



16(27r)^ 

+Ai(A;i, /c23)Ai(/c2, fci3)A3(A;3, A;i2) + Xi{ki, /c23)Ai(A;3, A;i2)A3(^2, ^13) 

(166) 



with 



A3(fci, ^23) = 7778^(34 - 5fc?)5(fc^ - A;23)e(fcf - k,) 

357r^ 

+-^(5(=')(k0(34 - 5kl,)QikF - hs) 

Kp 

= X3A{k,)^S{kF - ^23) + A3B(fc23)(2vr)^(5(ki) . (167) 
^23 

Integrating over the (5-functions leads to a lengthy expression that we will not 
give here. 

Using partial integration we can finally write the total expression in the form 

C[p] = qvpp - ^Cv2p . (168) 

The particular order of integrations we have carried out gives factors of kp 
appearing as UV cutoffs in the remaining integrals. Such a simplification arises 
for all contributions to the HF energy and the resulting integrals can there- 
fore be easily integrated numerically despite the relatively large number of 
integration variables. 

Key to the prescription used here is the Fourier transform of the expanded 
density matrices to momentum space. Due to its generality, this approach can 
easily be extended to the calculation of higher-order contributions to the DME. 
A similar approach was introduced in Ref. [15], where the authors used the 
Fourier transform of the expanded density matrix to generate medium inser- 
tions for a diagrammatic calculation of the nuclear energy density functional 
using chiral perturbation theory. 
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44-2 DME II 



The DME I prescription outlined above differs from tlie original NV approach 
in two respects. First, we do not rearrange and truncate the expansion by 
hand to ensure that the nuclear matter limit is exactly reproduced. Second, 
the DME 1 prescription keeps cross-terms in the product of the three expanded 
density matrices that arc formally of higher order in the NV approach. In order 
to quantify these effects and assess whether the expansion is under control, we 
have also performed the expansion where we strictly follow the original NV 
philosophy (DME II). 

We also note the differences in angle-averaging that arise with the different 
DME schemes. In the DME I approach, each p(xj, Xj) is first expanded in the 
natural Jacobi coordinates (R, rfe,ry), and then the three expanded density 
matrices are expressed in one common set of Jacobi coordinates. In the DME 
II prescription, we follow a different path by expressing the product of density 
matrices in one common set of Jacobi coordinates from the outset. The subse- 
quent DME implies a different angle-averaging, since only one density matrix 
is expanded in its natural Jacobi basis. We do not include the derivation of 
the DME II equations here, as it proceeds in much the same spirit as for the 
DME I, although we note that the final expressions are considerably more 
cumbersome since one finds different A; functions depending on whether one 
is expanding the p(xj, Xj) corresponding to the chosen Jacobi coordinates or 
one of the other two density matrices. 



5 Results 

In this section, we make some basic tests of the DME. We have two modest 
goals: to check that the DME does not degrade when applied to non-local, 
low-momentum NN potentials and to make a first assessment of the relative 
contributions of two- and three-body interactions. For the first goal, we ap- 
proximate the self-consistent Hartree-Fock ground-state wave function by a 
Slater determinant of harmonic oscillator single-particle wave functions. Us- 
ing these wave functions, we compare the DME approximation for the energy 
of a schematic model NN potential to the exact result where the finite range 
and non-locality of the interaction is treated without approximation. Then 
with the same wave function we check the error as we change the resolution 
(cutoff) of a realistic low-momentum potential. For the second goal, we ex- 
hibit some numerical results for the DME coefficient functions to illustrate 
the non-trivial density dependence and to show the effects of different pre- 
scriptions for the three-body DME. These are meant only to set a baseline 
because, at a minimum, we should include second-order contributions (i.e., 
beyond Hartree-Fock) before expecting quantitative predictions for nuclear 
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Fig. 4. Effects of different non-localities on the accuracy of the DME as a function 
of the ratio of the non-locality to range parameters for the harmonic oscillator 
approximation to the ground state of ^''Ca. 

structure or analyzing the cutoff dependence of the energy functional. How- 
ever, even at this stage it should be meaningful to use these results to compare 
the relative contributions of two- and three-body interactions. 

Although the original DME paper introduced formalism for non-local po- 
tentials [33], previous investigations of the effectiveness of the DME studied 
only local potentials (or local approximations to the G matrix). Because the 
low-momentum potentials used here can be strongly non-local, we first test 
whether the extra expansion required degrades the accuracy of the DME. We 
consider a model potential: 

with V a Gaussian potential, so the range is set by a. The range of the non- 
locality is set by /3; in the limit /3 — >• 0, V{y, r') — ^ v{v/a)5^{Y — r'). 

In Fig. m the effects of non-localities on the accuracy of the DME for in- 
tegrated quantities (e.g., {V)) is illustrated using this potential. We use a 
harmonic oscillator model of ^°Ca (i.e., the ground-state wave function is a 
Slater determinant of harmonic oscillator orbitals) and calculate the expec- 
tation value of the non-local V{y, r') in the Hartree-Fock ground state. For 
a given range a, we compare the error for a non-locality (3 to the error with 
/3 = 0. It is evident that the effect of the non-locality on the degradation of 
the DME is unimportant up to at least twice the range. Even when a is taken 
as small as the typical range of a repulsive core there should be no problem 
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Fig. 5. Errors per nucleon in the DME predictions for the expectation value of a 
model potential, Eq. (1169p . in a harmonic oscillator ground state for three N = Z 
nuclei (no Coulomb interaction). 

for the range of low-momentum cutoffs typically considered. 

The errors per nucleon for the DME with the same model ground state but 
with a realistic low-momentum nucleon-nucleon potential (starting from the 
chiral N^LO potential from Ref. ^9]) are shown in Fig. [5] for N = Z nuclei 
(without Coulomb) for A = 16, 40, and 80. It is evident that the cutoff de- 
pendence of the error is very slight until A < 2fm~^. Because the evolution 
of the potential does not alter the long-distance part, the weak cutoff depen- 
dence of the error implies that the short-distance contribution is very well 
reproduced and provides further confirmation that non-locality (which grows 
with decreasing A) is not a problem for the DME (note that long-range local 
interactions remain local). These errors are also smaller than errors found in 
early DME tests. 



The model calculations in Fig. [5] treat both direct (Hartree) and exchange 
terms with the DME. It was recognized long ago that the DME is ill-suited 
for long-range direct terms, which should be calculated exactly instead |34j . 
The dashed line in the figure shows the error for A = 40 but using the NLO po- 
tential, which does not have any long-range contributions to the direct scalar 
term. As expected, the error is significantly smaller than the N'^LO result, 
due at least in part to the crude treatment of the N^LO long-range direct 
contribution. Since the long-range local terms can be isolated in the poten- 
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Fig. 6. Contribution to the energy per particle in nuclear matter from the isoscalar 
coefficient function A{p) as a function of the density from the DME applied to the 
Hartree-Fock energy calculated using T^owfc with A = 2.1 fm~^. The result including 
the NN interaction alone is compared to NN plus NNN interactions for two DME 
expansions (I and II, see text). 

tial, it is feasible to perform exact Hartree evaluations of these pieces when 
implemented in a DFT solver. 

We turn now to the isoscalar A and B functions, which are the only con- 
tributors to uniform, symmetric nuclear matter. The energy per particle as a 
function of density p is given by: 



The individual contributions from A and B at the Hartree-Fock level are 
plotted in Figs. [6] and [71 and combined into E/A in Fig. [HI These use a two- 
body V\owk interaction evolved from the Argonne f is potential [70] with a sharp 
cutoff at A = 2.1 fm~^ and a chiral N^LO three-body force with constants fit 
to the binding energies of the triton and "^He [67] . Results are given using the 
NN contribution only and with NNN included, using the two prescriptions 
for the NNN double-exchange contribution (DME-I and DME-II) described in 
Section [H 

From Figs. [HI and [71 one sees that the ratios of contributions from three-body to 
two-body tend to increase monotonically with density, but are still only about 
20-30% at saturation density. This is consistent with general expectations 
from chiral power counting. The actual scaling with density of the ratio varies 
only slightly from being linear in the density. Because the local density in 
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Fig. 7. Contribution to the energy per particle in nuclear matter from the isoscalar 
coefficient function -B(p) as a function of the density from the DME applied to the 
Hartree-Fock energy calculated using Fiowfc with A = 2.1 fm~^. The result including 
the NN interaction alone is compared to NN plus NNN interactions for two DME 
expansions (1 and 11, see text). 
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Fig. 8. Energy per particle in nuclear matter by combining contributions from from 
the isoscalar coefficient functions Ai^p) and B(^p) with the kinetic energy as a function 
of the density from the DME applied to the Hartree-Fock energy calculated using 
Viowfc with A = 2.1 fm~^. The result including the NN interaction alone is compared 
to NN plus NNN interactions for two DME expansions (1 and 11, see text). Note 
that only expansion 11 correctly reproduces the nuclear matter limit. 
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actual nuclei in somewhat lower, there is reason to believe the expansion in 
many-body forces is under control. Past estimates of contributions to Skyrme 
energy functionals based on naive dimensional analysis [7Ij suggested large 
contributions from three-body and even four-body interactions. The present 
results imply more modest contributions, but evaluating the chiral N3L0 four- 
body contribution at Hartree-Fock will be needed for a definitive assessment. 

The comparison of the DME-I and DME-II curves gives us an estimate of the 
truncation error in the expansion applied to the NNN terms because these 
prescriptions differ in the contributions of higher-order terms in the expan- 
sion. Indeed, we have verified that suppressing these terms by hand brings 
the predictions for the A and B coefficients into agreement. The qualitative 
difference for the NNN-only contribution to B is large, but the actual coef- 
ficient itself is small, so this should not be alarming. However, because the 
combination of A and B and the kinetic energy to obtain the nuclear matter 
energy per particle involves strong cancellations, the spread in Fig. [8] is large 
on the scale of nuclear binding energies. 

These differences motivate a generalization of the Negele-Vautherin DME fol- 
lowing the discussion in Ref . |72j . In this approach, the expansion of the scalar 
density matrix takes the factorized form 

p(R + |, R - |) = ^ n„(A;ps) (a„(R)) , (171) 

n 

where 

(a„(R)) = {p(R), r(R), VV(R), ■ ■ ■ } , (172) 

and kp is a momentum scale typically taken to be kpCR,) as in Eq. fl52l) . Similar 
expansions are made for the other components of the density matrix. Input 
from finite nuclei can be used to determine the 11^ functions, which can be 
viewed as general resummations of the DME expansion; see Section [H] for a 
brief overview. 

Finally, in Fig. [9], the coefficient function C{p) is plotted as a function of 
density (p = 2A;p/37r^). Even at the highest density, the three-body contri- 
bution is a manageable correction to the two-body result. The NN -|- NNN 
result is in qualitative agreement with the results of Fritsch and collabora- 
tors who included two-pion exchanges with explicit A-isobars [17], although 
the three-body contributions in the current work are somewhat larger than ef- 
fects arising from explicit A-isobars. For this coefficient function the difference 
between DME-I and DME-II is comparatively small. 
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Fig. 9. Isoscalar coefficient function C{p) as a function of the density from the DME 
apphed to the Hartree-Fock energy calculated using 

Mowfc with A = 2.1 fm-^ The 
result including the NN interaction alone is compared to NN plus NNN interactions 
for two DME expansions (I and II, see text). 

6 Summary 

In this paper, we have formulated the density matrix expansion (DME) for 
low-momentum interactions and appUed it to a Hartree-Fock energy functional 
including both NN and NNN potentials. The output is a set of functions of 
density that can replace density-independent parameters in standard Skyrme 
Hartree-Fock energy density functionals. This replacement in Skyrme HF com- 
puter codes is shown schematically in Fig. (TDl Only one section of such a code 
would be replaced, and it takes the same inputs (single-particle eigenvalues 
and wave functions for the orbitals and the corresponding occupation num- 
bers) and delivers the same outputs (local Kohn-Sham potentials). Further- 
more, the upgrade from Skyrme energy functional to DME energy functional 
can be carried out in stages. For example, the spin-orbit part and pairing can 
be kept in Skyrme form with the rest given by the DME. Details of such a 
DME implementation will be given elsewhere. A further upgrade to orbital- 
based methods would also only modify the same part of the code, although 
the increased computational load will be significant. 

The numerical results given here are limited and do not touch on many of the 
most interesting aspects of microscopic DFT from low-momentum potentials. 
Topics to explore in the future include: 

• Examine the resolution or scale dependence of the energy functional by 
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evolving the input low-momentum potential. There will be dependence on 
the cutoff A (if using Vjowfc) or the flow parameter A (if using Vsrg) both 
from omitted physics and from intrinsic scale dependence. Calculations at 
least to second order are needed to separate these dependencies. 

• Examine the isovector part of the functional. We can isolate the contribu- 
tions from the more interesting long-range (pion) parts of the free-space 
interactions, allowing us to obtain analytic expressions for the dominant 
density dependence of the isovector DME coupling functions. 

• Study the dependence of spin-orbit contributions on NN vs. NNN inter- 
actions. This includes the isospin dependence as well as overall magni- 
tudes. The NN spin-orbit contributions arise from short-range interactions, 
whereas NNN contributions arise from the long-range two-pion exchange in- 
teraction. Therefore, we expect to find a rather different density dependence 
for the two types of spin-orbit contributions. 

• Explore the contribution of tensor contributions, which have recently been 
reconsidered phenomenologically [73|f7i] . 

• Understand the scaling of contributions from many-body forces. In partic- 
ular, how does the four-body force (which is known at N^LO in chiral EFT 
with conventional Weinberg counting) contribution at Hartree-Fock level 
impact the energy functional? 

The calculations presented here are only the first step on the road to a universal 
nuclear energy density functional (UNEDF) [15]. There are both refinements 
within the DME framework and generalizations that test its applicability and 
accuracy. While many of these steps offer significant challenges, in every case 
a plan is in hand to carry it out. The DME can be directly extended to 
include second-order (or full particle-particle ladder) contributions by using 
averaged energies for the energy denominators. However, a more systematic 
approximation is under development using a short-time expansion [75j. More 
difficult future steps include dealing with symmetry breaking and restoration 
in DFT for self-bound systems, dealing with non-localities from near-on-shell 
particle- hole excitations (vibrations), and incorporating pairing in the same 
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energy 
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Fig. 10. Diagrams showing the flow in Skyrme HF codes at present (left) and mod- 
ified for the DME (right). 
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microscopic framework (see Ref. [H]). 



In extending our calculations we will also modify the standard DME formalism 
from Ref. [33j that we have followed in the present work. The formalism has 
problems even beyond the truncation errors from different DME prescriptions 
already discussed in Sections H] and O the most severe being that it provides an 
extremely poor description of the vector part of the density matrix. While the 
standard DME is better at reproducing the scalar density matrices, even here 
the errors are sufficiently large that the disagreement with a full finite-range 
Hartree-Fock calculations can reach the MeV per particle level. Gebremariam 
and collaborators have traced both of these problems to an inadequate phase 
space averaging (PSA) used in the previous DME approaches [39]. In the 
derivation of the DME, one incorporates average information about the local 
momentum distribution into the approximation. The Negele-Vautherin DME 
uses the phase space of infinite nuclear matter to perform this averaging. 
However, the local momentum distribution in finite Fermi systems exhibits 
two striking differences from that of infinite homogenous matter. First, mean- 
field calculations of nuclei show that the local momentum distribution exhibits 
a diffuse Fermi surface that is especially pronounced in the nuclear surface. 
Second, the local momentum distribution is found to be anisotropic, with the 
deformation accentuated in the surface region of the finite Fermi system. 

To incorporate both of these missing effects into the DME, Gebremariam et 
al. have constructed a model for the local momentum distribution based on 
previous studies of the Wigner distribution function in nuclei [2H]- The model 
parameters are adjusted so that the DME accurately reproduces both inte- 
grated quantities, such as the expectation value of the finite-range nucleon- 
nucleon interaction taken between Slater determinants from self-consistent 
Skyrme-Hartree-Fock calculations, as well as the density matrices themselves. 
The improvements are substantial, typically reducing relative errors in inte- 
grated quantities by as much as an order of magnitude across many different 
isotope chains. The improvement is especially striking for the vector density 
matrices. We will test this improved DME in future investigations. 

The tests of the DME will include benchmarks against ab initio methods in 
the overlap region of light-to-medium nuclei. Additional information is ob- 
tained from putting the nuclei in external fields, which can be added directly 
to the DFT/DME functional. Work is in progress on comparisons to both 
coupled cluster and full configuration interaction calculations. A key feature 
is that we use the same Hamiltonian for the microscopic calculation and the 
DME approximation to the DFT. The freedom to adjust (or turn off) external 
fields as well as to vary other parameters in the Hamiltonian permits detailed 
evaluations of the approximate functionals. In parallel there will be refined 
nuclear matter calculations; power counting arguments from re-examining the 
Brueckner-Bethe-Goldstone approach in light of low-momentum potentials 
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will provide a framework for organizing higher-order contributions. These in- 
vestigations should provide insight into how the energy density functional can 
be fine tuned for greater accuracy in a manner consistent with power counting 
and EFT principles. 
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